function out = runOptimization(setupFilter,params)
%% For logging
if setupFilter.MEXon == 1
    if setupFilter.optim > 0 && setupFilter.numCPUs > 2
        now             = clock;
        nameOfLogFile   = ['logFile_month',num2str(now(2)),'_day',num2str(now(3))...
            ,'_hour',num2str(now(4)),'_min',num2str(now(5)),'.txt'];
        %diary(nameOfLogFile)
    end
end
%% The estimation
paramsValues           = struc2values(params,setupFilter.selectParams);
exitflag               = 0;
if setupFilter.optim == 1
    options = optimset('Display','iter','MaxIter',setupFilter.MaxIter,...
        'MaxFunEvals',setupFilter.MaxEvals,'TolFun',setupFilter.TolFun,'TolX',setupFilter.TolX);
    [paramsOptValues,fval,exitflag] = fminsearch(@objectFuncQML_ForOptim,paramsValues,options,setupFilter);    
elseif setupFilter.optim == 2
        % For Nelder-Mead with transformed space for the parameters
    options = optimset('Display','iter','MaxIter',setupFilter.MaxIter,...
        'MaxFunEvals',setupFilter.MaxEvals,'TolFun',setupFilter.TolFun,'TolX',setupFilter.TolX);
    setupFilter.transformParamsOn = 1;
    paramsValues    = transformParams(paramsValues,setupFilter.lowerBoundsValues,setupFilter.upperBoundsValues);
    [paramsOptValues,fval,exitflag] = fminsearch(@objectFuncQML_ForOptim,paramsValues,options,setupFilter);
    paramsOptValues = unTransformParams(paramsOptValues,setupFilter.lowerBoundsValues,setupFilter.upperBoundsValues);    
    setupFilter     = rmfield(setupFilter,'transformParamsOn');   
elseif setupFilter.optim == 3
    % Matlab's gradient-based optimizer: Does not work well!
    %options = optimset('Display','iter','MaxIter',setupFilter.MaxIter,...
    %    'MaxFunEvals',setupFilter.MaxEvals,'TolFun',setupFilter.TolFun,'TolX',setupFilter.TolX,...
    %    'GradObj','on','Hessian','user-supplied','Algorithm','trust-region-reflective');
    %[paramsOptValues,fval,exitflag,output,lambda,grad,hessian]= fmincon(@objectFuncQML_ForOptim,paramsValues,...
    %    [],[],[],[],setupFilter.lowerBoundsValues,setupFilter.upperBoundsValues,[],...
    %    options,setupFilter);
    % Own implementation of the Newton-optimizer
    setupFilter.transformParamsOn = 1;
    paramsOptValues = newtonOptim(paramsValues,setupFilter,...
        setupFilter.MaxIter,setupFilter.TolFun,setupFilter.TolX);   
    setupFilter = rmfield(setupFilter,'transformParamsOn');

elseif setupFilter.optim == 4
    % the CMAES
    InsigmaValues      = setupFilter.InsigmaValues;          %The std deviation in the initial search distributions
    sigma              = setupFilter.cmaesSigma;             %The step size
    opts.SigmaMax      = setupFilter.cmaesSigmaMax;          %The maximal value for sigma
    opts.LBounds       = setupFilter.lowerBoundsValues;      %Lower bound for params
    opts.UBounds       = setupFilter.upperBoundsValues;      %Upper bound for params
    opts.MaxIter       = setupFilter.MaxIter*100;             %The maximum number of iterations
    opts.MaxFunEvals   = setupFilter.MaxEvals*100;            %The maximum number of function evaluations
    opts.PopSize       = setupFilter.cmaesPopSize;           %The population size
    opts.VerboseModulo = 1;                                  %Display results after every 10'th iteration
    opts.TolFun        = setupFilter.TolFun;                 %Function tolerance
    opts.TolX          = setupFilter.TolX;                   %Tolerance in the parameters
    opts.Plotting      = 'off';                              %Dislpay plotting or not
    opts.Saving        = 'off';                              %Saving results
    opts.numCPUs       = setupFilter.numCPUs;
    %paramsOptValuesCMEAS = cmaes_dsgeDisplay_mp(@objectFuncQML_ForOptim,paramsValues,sigma,InsigmaValues,opts,setupFilter);    
    %paramsOptValues = paramsOptValuesCMEAS;
    [xmin,fmin,counteval,stopflag,outhist,bestever] = cmaes_dsgeDisplay_mp_grendel(@objectFuncQML_ForOptim,paramsValues,sigma,InsigmaValues,opts,setupFilter);    
    paramsOptValues = bestever.x;
elseif setupFilter.optim == 5
    options = optimset('Display','iter','MaxIter',setupFilter.MaxIter,...
        'MaxFunEvals',setupFilter.MaxEvals,'TolFun',setupFilter.TolFun,'TolX',setupFilter.TolX);
    [paramsOptValues] = fmincon(@objectFuncQML_ForOptim,paramsValues,[],[],[],[],...
        setupFilter.lowerBoundsValues,setupFilter.upperBoundsValues,[],options,setupFilter);
else
    paramsOptValues = paramsValues;
end
paramsOpt = values2struct(paramsOptValues,setupFilter.selectParams);
[Q,errorMes,model,setupFilter,outFilter] = objectFuncQML(paramsOptValues,setupFilter);
if errorMes == 1
    out = [];
    return
end

% We always want to have the model fit for surveys computed
modelPosterior = model;
modelPosterior.inclSurveys = 1;
if model.inclSurveys == 0
   % Solve for expected short rates
   var0   = model.by0(1,1);
   varx   = model.byx(1,:);
   varxx  = model.byxx(1,:,:);
   varxxx = model.byxxx(1,:,:,:);
   [r0,rx,rxx,rxxx]= projectionStepExpectedControlClosedForm(var0,varx,varxx,varxxx,...
       max(setupFilter.matConShortRate),model,setupFilter.setupProj);
   modelPosterior.matConShortRate = setupFilter.matConShortRate;
   modelPosterior.r0     = r0;
   modelPosterior.rx     = rx;
   modelPosterior.Rxx    = reshape(rxx,length(setupFilter.matConShortRate),model.nx^2);
   modelPosterior.Rxxx   = reshape(rxxx,length(setupFilter.matConShortRate),model.nx^3);
end   
version        = 1+1;   %1 for evaluate at xhat and no account of state uncertainty
                        %2 for evaluate at xhat and accounting for uncertainty 
dimy           = setupFilter.numMacro+length(setupFilter.maturities)+length(setupFilter.matConShortRate);                         
outFilter.yhat = CDKF_y_posterior(setupFilter.data(:,1:dimy)',modelPosterior,setupFilter,outFilter.xhat,outFilter.Smat,version);
outFilter.what = outFilter.yhat - setupFilter.data(:,1:dimy)';

% Extra yields
outFilter.YieldsAll_cu = 4*getYields(modelPosterior,outFilter.xhat);

% Add a simulated sample path to the reported output if not already present
if setupFilter.numSim == setupFilter.burnIn
    setupFilterSim = setupFilter;
    setupFilterSim.numSim = size(setupFilterSim.shocks,2);
    [QSim,~,~,setupFilterSim,outFilterSim] = objectFuncQML(paramsOptValues,setupFilterSim);
    outFilter.outSim                = outFilterSim.outSim;
    outFilter.modelMomentsSim       = outFilterSim.modelMoments;
    outFilter.xSim                  = outFilterSim.outSim.x_cu;
    outFilter.yieldsAllSim          = outFilterSim.yieldsAllSim;
    outFilter.yieldsAllSimWithNoise = outFilterSim.yieldsAllSimWithNoise;
    outFilter.ySimAsInData          = outFilterSim.ySimAsInData;
    outFilter.ySimAsInDataWithNoise = outFilterSim.ySimAsInDataWithNoise;
    outFilter.QSim                  = QSim;
end

% Adding term premia for nominal yields to model
model        = getTermPremia_loadings(model,setupFilter);
termPrem     = getTermPremia(model,outFilter.xhat);
[termPremSim,TPsim]  = getTermPremia(model,outFilter.outSim.xAll_cu);
T            = size(termPrem.TP,2);
TPdecom      = zeros(length(model.by0),T,model.ne);
modelTmp     = model;
for i=1:model.ne
    % Version 1
    %xhatDecom = 0*outFilter.xhat;
    %xhatDecom(1:model.nx1,:) = outFilter.xhat(1:model.nx1,:);   %endogenous states
    %xhatDecom(model.nx1+i,:) = outFilter.xhat(model.nx1+i,:);   
    %tmp = getTermPremia(modelTmp,xhatDecom);
    %TPdecom(:,:,i) = tmp.TP;
    
    %Version 2 (should be more accurate than version 1)
    tmpAll = getTermPremia(modelTmp,outFilter.xhat);
    xhatDecom = outFilter.xhat;
    xhatDecom(model.nx1+i,:) = 0*outFilter.xhat(model.nx1+i,:);   %turn-off shock i
    tmp = getTermPremia(modelTmp,xhatDecom);
    TPdecom(:,:,i) = tmpAll.TP-tmp.TP;
end
termPrem.TPdecom = TPdecom;

% Computing real yields
setupFilter.setupProj.objectFuncNum = 3; %for REAL bond prices
modelRealYields      = model;
modelRealYields      = rmfield(modelRealYields,'thetaBondPrice');
modelRealYields      = projectionStepBondPricesClosedForm(modelRealYields,setupFilter.setupProj);
realYields           = 4*getYields(modelRealYields,outFilter.xhat);
realYieldsSim        = 4*getYields(modelRealYields,outFilter.outSim.xAll_cu);
% Adding term premia for real yields to modelRealYields
modelRealYields      = getTermPremia_loadings(modelRealYields,setupFilter);
[termPremRealSim,TPrealSim] = getTermPremia(modelRealYields,outFilter.outSim.xAll_cu);
termPremReal         = getTermPremia(modelRealYields,outFilter.xhat);
inflPremSim.meanTP   = mean(TPsim-TPrealSim,2)*40000;
inflPremSim.stdTP    = std(TPsim-TPrealSim,0,2)*40000;
% The equity price
setupFilter.setupProj.objectFuncNum = 5; %for Equity price
consClaim      = 1;
model          = projectionStepEquityPrice(model,setupFilter.setupProj,consClaim);
consClaim      = 0;
model          = projectionStepEquityPrice(model,setupFilter.setupProj,consClaim);
outSim         = simProjection(model,setupFilter.shocks);
outHat         = getAllControls(model,outFilter.xhat);
consClaim      = 1;
equityPremConsSim = getEquityPrem(model,modelRealYields,outSim,setupFilter.T,consClaim);
equityPremConsHat = getEquityPrem(model,modelRealYields,outHat,setupFilter.T,consClaim);

consClaim         = 0;
equityPremSim     = getEquityPrem(model,modelRealYields,outSim,setupFilter.T,consClaim);
equityPremHat     = getEquityPrem(model,modelRealYields,outHat,setupFilter.T,consClaim);

equityPremData    = getEquityPremData(setupFilter);

% Risk-adjusted CS regression loadings
% Convert into forwards
fwdPrem = nan(size(termPrem.TP));
for i = 1:size(termPrem.TP,1)-1
     fwdPrem(i,:) = (i+1)*termPrem.TP(i+1,:) - i*termPrem.TP(i,:);
end
% TPxhat and hence fwdPrem are scaled by 40000 to get risk premia in
% annualized basis points in getTermPremia. 
% Annualized yields in are not scaled. Hence, we need to undo the scaling of TP 
% to get basis points, but we still need to preserve the annualization i.e. the multiplication
% by 4.
numBoots = 10000;
CS_riskAdj = CampbellSchillerRegression_riskAdjust(setupFilter.yieldsAll,...
     termPrem.TP'/40000*4,fwdPrem'/40000*4,setupFilter.CSregress_m,numBoots);
%% Output and standard errors
out.exitflag          = exitflag;
out.paramsOpt         = paramsOpt;
out.model             = model;
out.modelRealYields   = modelRealYields;
out.realYields        = realYields;
out.realYieldsSim     = realYieldsSim;
out.setupFilter       = setupFilter;
out.outFilter         = outFilter;
out.termPrem          = termPrem;
out.termPremSim       = termPremSim;
out.termPremReal      = termPremReal;
out.termPremRealSim   = termPremRealSim;
out.inflPremSim       = inflPremSim;
out.equityPremConsSim = equityPremConsSim;
out.equityPremConsHat = equityPremConsHat;
out.equityPremSim     = equityPremSim;
out.equityPremHat     = equityPremHat;
out.equityPremData    = equityPremData;
out.CS_riskAdj        = CS_riskAdj;
epsValue = 1D-5;
% Compute standard errors if possible
if any(paramsOptValues-epsValue< setupFilter.lowerBoundsValues) ...
	|| any(paramsOptValues+epsValue > setupFilter.upperBoundsValues)
    out.SEcanBeComputed = 0;
    return;
else
    out.SEcanBeComputed = 1;
    setupFilter.setupProj.dispOn = 0;
    if setupFilter.seOn == 1
        for i=5:7
            setupFilter.epsValue = 10.^(-i);
            tmp = num2str(setupFilter.epsValue);
            if setupFilter.epsValue == 1D-2
               name = 'paramsSE_eps2';
            elseif setupFilter.epsValue == 1D-3
               name = 'paramsSE_eps3';
            elseif setupFilter.epsValue == 1D-4
               name = 'paramsSE_eps4';
            elseif setupFilter.epsValue == 1D-5
                name = 'paramsSE_eps4';
            else
               name = ['paramsSE_eps',tmp(5)];
            end
            disp([name,['  epsValue = ',num2str(setupFilter.epsValue)]])
            [out.(name)] = getQMLse(paramsOptValues,setupFilter,outFilter);
        end
    end
end
if setupFilter.MEXon == 1
    if setupFilter.optim > 0
        %diary off
    end
end

end